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Data driven optimal filtering for phase and frequency of noisy oscillations: application 

to vortex flowmetering 
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A new method for extracting the phase of oscillations from noisy time series is proposed. To 
obtain the phase, the signal is filtered in such a way that the filter output has minimal relative 
variation in the amplitude (MIRVA) over all filters with complex-valued impulse response. The 
argument of the filter output yields the phase. Implementation of the algorithm and interpretation 
of the result are discussed. We argue that the phase obtained by the proposed method has a low 
susceptibility to measurement noise and a low rate of artificial phase slips. The method is applied 
for the detection and classification of mode locking in vortex flowmeters. A novel measure for the 
strength of mode locking is proposed. 

PACS numbers: 05.45.Tp, 06.30.Ft, 05.45.Xt 
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I. INTRODUCTION 

Several modern methods for time series analysis make 
explicit use of the phase of measured oscillatory signals. 
Examples are tests for unidirectional [l| or mutual 
synchronization of chaotic oscillators, based on accurate 
or noisy data, identification of the coupling direction 
0,0, or indicators for generalized synchronization 0. 
Phase analysis has successfully been applied in neurology 
P, cardiology HI IH Ell i ecology El > an d astronomy 
|12| (for recent, comprehensive reviews, see 0,0]). 

A number of methods have been proposed for the de- 
termination of the phase. Among these are (a) phase 
extraction from the argument of the analytic signal |15| , 
from the convolution of the signal with a Morlet wavelet 
|l6l IT?] ], or after complex demodulation or quadrature 
filtering 0, (b) the angle of circulation of a 2D projec- 
tion of the reconstructed phase-space trajectory 0| or 
its time derivative [l9l | around a point, and (c) linear in- 
terpolation 1] of phase between distinct events marking 
the beginning new cycles [l4| . 

Although some rules for selecting the appropriate 
method for a given system have been proposed the 
choice is not always obvious. The wish list for proper- 
ties of the reconstructed phase 4>(t) includes: a con- 
stant advance of 2w per cycle, a steady accumulation 
{(j>(t) « const.), accuracy in the presence of measure- 
ment noise, unambiguity with respect to 2n phase slips, 
and a functional dependence on the current state of the 
oscillator (locality). Autonomy of the oscillator combined 
with steady accumulation and locality of phase implies 
that 4>(t) is the variable that corresponds to the zero Lya- 
punov exponent of the system - another desired property. 

But only for perfectly periodic signals can all these 
wishes be fulfilled. For deterministic, chaotic oscillators 
the linear-interpolation method (c) does often lead to 
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a satisfactory steady and local phase. The problem of 
defining a steady, local phase when the internal dynam- 
ics of a deterministic, chaotic oscillator are known was 
treated rigorously in Ref . j2£j • But as internal- and mea- 
surement noise become stronger, some temporal averag- 
ing is required and locality in time has to traded for accu- 
racy and/or unambiguity. In Refs. 0, IS H> El El 
the condition of unambiguity was relaxed and only the 
cyclic phases [4>{t) mod 27r] of noisy oscillations were used. 
Thus data analysis was insensitive to phase slips, i.e., 
sudden advances of the phase by ±27r, which may or may 
not be artifacts of measurement noise. Here we choose to 
be less demanding with respect to locality, in favor of a 
steady, accurate and, as much as possible, unambiguous 
phase. 

In order to identify a corresponding method of phase 
extraction, notice that the computation of the analytic 
signal (or complex demodulation) [method (a)] is gener- 
ally recommended to be combined with linear band-pass 
filtering of the desired oscillatory component 0, . 
The overall effect is the application of a complex- valued, 
linear band-pass filter When the method of delays 
is used for the phase-space reconstruction of the angle-of- 
circulation method (b), the 2D projection is also equiv- 
alent to complex- valued linear filtering |22j . likewise for 
calculation of time derivatives. 

Finally, in the vicinity of a Hopf bifurcation, where 
dynamics can be brought into Hopf normal form by a 
nonlinear coordinate transformation (see, e.g., [23] i. this 
transformation is done in such a way that all contribu- 
tions to dynamics which are "non-resonant" with the os- 
cillation at the fundamental frequency are eliminated. In 
mere kinetic terms this simply amounts to eliminating 
higher harmonics and offsets, which can be achieved by 
complex, linear band-pass filtering. The SU(2) symme- 
try of the Hopf normal-form guarantees the steady ac- 
cumulation of phase in the steady state, when phase is 
measured as the angle of rotation around the origin. 

Thus, a unified view on (a) and (b) is complex, linear 
band-pass filtering. It achieves steady accumulation in 
a natural way when the fundamental mode is isolated. 



The results regarding accuracy and unambiguity depend 
on the choice of the filter. Since the concept of phase orig- 
inates from limit-cycle oscillations, which, in the trans- 
formed coordinates, correspond to a motion on a circle, 
our idea for choosing the filter here is to make the filtered 
signal move as close as possible to a circle in the com- 
plex plane. Roughly speaking, we consider the motion on 
the circle as the signal and the deviations as noise and 
maximize the signal to noise ratio (SNR) - even though 
not all deviations are actually due to measurement noise. 
Since, with such a filter, noise-induced excursions of the 
trajectory to the origin of the complex plane are mini- 
mized, this is also a good way to reduce ambiguities in 
the phase. The maximization of the SNR is done not 
only with respect to width and center frequency of the 
filter, but with respect to the complete dynamics of its 
impulse response. The determination of the filter is non- 
parametric and data driven. 

We proceed as follows. Section [D] contains a mathe- 
matical formulation of the ideas outlined above and lists 
some implications. In Section IIIII the method is applied 
to simulated data, with special attention to the effect of 
filtering on measured frequencies. A practical applica- 
tion to vortex flowmeters is discussed in Sec. IIVI where 
we also introduce a novel measure for the strength of 
mode locking. 



A(0=8.0 



A(0=1.0 



MIRVA 



1 

-1 



b) 

9 

s 6 
a> 3 




V - 





w original 
Aq>=8.0 
Aq>=1.0 
MIRVA 



100 



200 



300 



FIG. 1: Three different filters applied to the same time 
series 10 • a) The demodulated, filtered signals Z(t) = 
z(t) exp(— iu>o t). b) The phase <j>o(t) — ^ot of the original 
signal and the relative phases (f>(t) — 4>o(t) obtained using the 
three filters. See text for details. 



2. Example 



II. THEORY 



MIRVA niters 



As an example, assume that x(t) is composed of a 
constant-amplitude oscillation with phase fluctuations 
and white measurement noise, 



x{t) = cos[w t + <j) Q (t)} + l](t), 



(2) 



1. Definition 

Let x(t) be a real- or complex- valued stationary signal 
with oscillatory components. Denote by z(t) the signal 
obtained from x(t) by linear filtering with a complex- 
valued filter with impulse response f(t), i.e., z = f * x. 
Define q as the nonnegative number such that 
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where (•) denotes the expectation value. 

Now, search for a filter / such that the quantity q 
given by Q has a local minimum with respect to the 
filter. Such a filter / minimizes the relative variance of 
the amplitude (MIRVA) for the given signal x(t). The 
practical computation of MIRVA filters is addressed in 
Appendix El 

For every MIRVA filter /, there is a two-parameter 
family of MIRVA filters f s ,c{t) = cf(t — s) with real s and 
complex c. Below we shall always have a single member 
stand for the whole family without saying. 



where 



(Mt)Mt')) = 2DS(t - 0, WW)) = 2G5(t - 



0- 

(3) 



This signal mimics narrow-banded limit-cycle and 
chaotic oscillations in the vicinity of the fundamental 
frequency. The larger the noise strength G, the more 
difficult the determination of the phase <fio(t) from x(t) 
becomes. We simulate x(t) with D = 1, G = 0.01, and 
ujq = 20 over an interval of length T — 256. Without any 
filtering, the SNR is zero. Figure^, shows the demod- 
ulated signals Z(t) — z[t) exp(— iu>o t) for three different 
complex filters fit)- The first two filters are of the form 



f(t) - exp 



(4) 



One is a comparatively wide (Aw = 8.0) band pass, the 
other one is rather narrow (Aw = 1.0). The third is the 
MIRVA filter obtained by minimizing Q. It is approx- 
imated by @ with Auj — 2.9. As is shown in Fig. ^d, 
both the narrow and the wide filter lead to artificial phase 
slips. Only when using the MIRVA filter is 4>o(t) faith- 
fully tracked. 
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3. Remarks on the minimization of q 

Let q m i n denote the value of q attained at a local mini- 
mum. Since the operation of linear filtering defines a semi 
group, g m i n is also a local minimum of q with respect to 
further filtering of z(t), i.e., for q calculated with z re- 
placed by z' := g * z. The minimum is attained when the 
filter g is the unit element of the semi group, the Dirac 
(5-function. As a result one has 
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for all r. In particular, when differentiating with respect 
to r at r = and taking the imaginary part, it follows 
that 
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where u>i is the instantaneous frequency defined by 
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and Wmoan is known as the mean frequency, defined either 
by the last equation of © or, equivalently, as the "center 
of mass" of the power spectrum of the signal z. For the 
relation of the mean frequency to the phase frequency (or 
average frequency) w p h := (u>i) see Sec. IIIII 

Often, signals contain oscillations at several different 
frequencies. A systematic method to extract various os- 
cillation frequencies has been proposed in [24| . When us- 
ing the concept of MIRVA filtering, distinct local minima 
of q can be identified with distinct oscillatory components 
of the signal. 



4- Special cases 

In two special cases the problem of finding MIRVA 
filters can be discussed analytically: For perfectly peri- 
odic or quasi-periodic signals there is, for every Fourier 
mode excited by the signal, a MIRVA filters that ex- 
tracts exactly this mode. The filtered signal is of the 
form z(t) = expiicot) and q m in — 0. This holds true also 
if the signal is overlaid by any kind of noise. 

For Gaussian, linear processes it is always possible to 
find filters such that (z 2 ) = and (\z\ 2 ) ^ 0, e.g., by 
letting only Fourier modes with positive frequency pass. 
Then (|z| 4 } = 2 (|z| 2 ) 2 , q = q m - m = 1, and all these filters 
are MIRVA filters. 

The processes we are interested here are typically lo- 
cated between these two poles: Noisy, nonlinear, periodic 
processes with some fluctuations in the phase. Thus we 



expect < <7 m i n < 1. When solving the optimization 
problem q 2 = min numerically with time series of finite 
length T (see Appendix |A"|> . local minima with g m j n > 1 
have also been found. It is not clear whether these persist 
in the limit T —> oo. 



B. The phase of MIRVA filtered signals 

1. Definition and error estimates 

The main purpose of MIRVA filtering is to obtain the 
phase 

4>{t) = J 0Ji(t')dt' = argz (mod27r) (8) 

of the oscillations extracted by the filter. 

Denote by 77(f) the contribution of measurement noise 
to z(t), the "true" value of z(t) by zq (f) := z(t) — 77(f), 

and the "true" phase by (f>o(t) := / Im{zo(t') / Zo(t')}dt' . 
Two kinds of errors in (j)(t) caused by measurement noise 
can be distinguished: deviations by multiples of 2n, i.e. 
phase slips, which are due to noise-induced excursions of 
z(t) around the origin and accumulate as time proceeds, 
and errors in the cyclic phase [0(f) — (f>o(t) + ir] mod27r, 
which have a finite correlation time. The distinction is 
particularly sharp when q min is small enough, so that the 
probability density for values of z(t) near zero is small, 
or, as we shall consider now, when q is small for general 
filters /. 

An order of magnitude estimate for an upper bound to 
the rate of noise-induced phase slips is given by 



P 



<M 2 > 



= 



■Au 



(9) 



where the first term denotes the probability density of 
|z| 2 /(|z| 2 ) at zero and the second term denotes the spec- 
tral width of the filter. The first term typically decays 
exponentially fast as q 2 decreases, while the relation be- 
tween q 2 and Aw is only algebraic in general. Hence, 
minimizing q 2 is a good strategy to minimize phase slips. 

For the noise-induced error in the cyclic phase, an ex- 
act upper bound can be obtained in the limit that q 2 is 
small. For simplicity, assume that the noise has under- 
gone sufficient temporal averaging by the filter /, so that 
the central limit theorem applies and 77(f) is Gaussian. 
Since / is a complex band-pass filter, (?7 2 ) = 0. In order 
to derive an upper bound for (I77I 2 ) from q 2 , we assume 
the worst case, that is, all variation in \z\ 2 is due to 77 
only, while |zo| 2 = const. = 1 without loss of generality. 
Since 77 is independent of zq, Eq. Q can then be written 
as 



2\2 



1 + 4(|77| 2 )+2(|77| 2 ) 
1 + 2<M 2 ) + (|77| 2 } 2 



(10) 
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Solving for (\f]\ 2 ) yields 



The corresponding noise-induced variance in <f> is 
var^> = var [arg(z +ry) — arg(z )] 



; var 



; var 



arg ( 1 + — 



+ 0((\r 1 \ 2 ) 2 ) 



(\v(t)\ 2 ) 



2\2x 



o((\v\ 2 ) 



(11) 



(12) 



In the general case, when \zo\ is not constant, we get, 
from combining pif) and (|12[) . 



var d> < 



0(q 4 ). 



(13) 



in a sufficiently high-dimensional delay space through a 
particular kind of Poincare section called a counter. For 
the topological frequency, the invariance with respect to 
linear filtering has been proven rigorously [26j. When 
the modulus of the signal Zi obtained by MIRVA filtering 
a signal Xi and its linear interpolation have a nontrivial 
lower bound, i.e., \l Zj + (1 — 1) > d > for < / < 1, 
and when the impulse response fk of the MIRVA filter de- 
cays sufficiently fast for large k, then w p h obtained from 
Zi and its linear interpolation is (up to the sign) identi- 
cal to a topological frequency of the oscillations of the 
signal Xi. A corresponding counter can be obtained as 
follows: Assume that all significant contributions to f^ 
are within a range of M consecutive delay times. Then 
the filter operation f * x can be interpreted as a projec- 
tion from the M-dimensional delay space of Xi into the 
two-dimensional complex plane. The counter is given by 
all points in delay space which are projected onto the 
real, non-negative half axis. 



Thus, minimizing q 2 is a good strategy to minimize noise- 
induced errors in the measured cyclic phases. 



III. THE EFFECT OF MIRVA FILTERING ON 
MEASURED FREQUENCIES 



2. Phase diffusion 

Over long time intervals, <p(t) typically performs a ran- 
dom walk with drift. Thus, another important character- 
istic of the phase is its diffusion coefficient 

D := hm W + T)- m-^T) 

T — >oo 2 1 

The estimation of D from finite-length samples of <p(t) is 
discussed in Appendix IbI 



C. Invariance with respect to filtering of x(t) 

The MIRVA filtered signal z(t) and the phase and fre- 
quency derived thereof are invariant with respect to lin- 
ear filtering of the original signal x(t) in the following 
sense: Let y{t) be a signal obtained from x(t) by linear 
filtering, i.e. y = h * x, and let / be a MIRVA filter for 
x. Then, a MIRVA filter for y is, at least formally, given 
by /' = / * h , where h^ 1 is the inverse of h defined 
by hr 1 * h = 5. So the filtered signal z = f*x = f'*y 
which satisfies the minimization condition (J5J is identical 
for x and y. For example, MIRVA filtering can be used to 
extract the phase and frequency of an oscillatory signal 
which has been "bleached" |25| . i.e., filtered such as to 
make its power spectrum white (see |26|L 

The concept of MIRVA filtering carries straightfor- 
wardly over to a discrete-time representation of signals 
Xi, Zi and filters sampled at time intervals of length 
At. In Ref. j2|| the notion of a topological frequency of a 
time-discrete signal xi is defined. Roughly speaking, this 
is the rate of transitions of the trajectory of the signal 



In order to illustrate the effects of MIRVA filtering on 
measured frequencies, the method is applied to a numer- 
ical solution of the noisy Stuart-Landau Equation (or 
Hopf normal form) 



A = (e + iuji)A -(g r + i gi)\A\ 2 A + £ 



(15) 



where A = A(t) represents the complex amplitude of 
an oscillator and £(£) is complex- valued, Gaussian, white 
noise with correlations 

(C(«X(*')> = and (mat')) = 4GS(t - t'). (16) 

In a certain sense, this system universally describes noisy 
oscillations in the vicinity of a Hopf bifurcation |2jj ■ As- 
sume the bifurcation to be supercritical (g r > 0) and do 
a linear change of coordinates to set g r = G = 1 and 
lui = without loss of generality (even though, in prac- 
tice, uji ^> e). When gi ^ 0, the phase frequency, 



j pKA — —<)i 



■2Af- 



(17) 



(where N := tt 1 / 2 exp(e 2 /4)[l + erf(e/2)], see, e.g., jH 
|28|l calculated from lo^a = lra{A/A} directly without 
filtering, differs from the corresponding mean frequency 



-9i 



2C 



e-4 2e 



2 My 1 



(18) 



This is a direct consequence of the correlation between 
uh,a and |^1| 2 ((u; hA \A\ 2 ) / (\A\ 2 ) ? (u> i>A )). 

As an example, a simulation of A(t) of length T = 10 6 
with e = 2, and gi = 1 is generated using Euler steps 
of length 2~ 11 . For these parameters, ui p h,A — —2.225, 
w me an,A = —2.899, and the relative variance of the unfil- 
tered amplitude is (|A| 4 ) / (\A\ 2 ) 2 - 1 = (2A/" 2 - 2eAf - 
4)/(eA/"+ 2) 2 = 0.303. 
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FIG. 2: MIRVA filtering of the process A(i) given by 
Eqs. I|15llf)|l with lui = 0, g r = g% = G = 1, and e = 2. 
The lower graph is a blow-up of the upper graph. Both show 
the power spectral density of A(i) (solid), of the filtered sig- 
nal z(t) (dotted), and the characteristics of the MIRVA filter 
(dash-dotted). The vertical lines indicate the locations of var- 
ious frequencies associated with the process. 



The MIRVA filter for the time series x{t) = A{t) is cal- 
culated by the indirect method described in Appendix lAl 
The filter reduces the relative variance of the ampli- 
tude to g^ in = 0.164. Figure [21 shows the characteris- 
tics of the MIRVA filter in comparison with the power 
spectrum of the original signal A(t) and the filtered sig- 
nal z(t). The locations of phase- and mean frequency 
before filtering and after filtering [w p h, z = —2.245(2), 
Wmean,z = -2.298(2)] are also indicated. The MIRVA 
filter is a rather wide (half-width 2Au « 6.5), approx- 
imately symmetric band-pass filter with a center fre- 
quency lu c w —1.3 below the linear frequency u>i = 0. 
As a result, the phase frequency of the filtered sig- 
nal is also shifted to lower frequencies, but the effect 
w P h,z - w P h,A ~ -0.02 = 0(Auj- 1 Af~ 1 LJ c ) (see 0) is 
quite small. 

On the other hand, there is a pronounced shift in the 
mean frequency by MIRVA filtering: the mean frequency 
approaches the phase frequency. This is a generic effect 
of MIRVA filtering. In the presence of a negative corre- 
lation between amplitude and instantaneous frequency, 
as found in our example, a filter that amplifies the sig- 
nal when LOi is high and damps the signal when Ui is low 
reduces fluctuations in the amplitude and, at the same 
time, reduces the correlation. Since the phase frequency 
is only weakly affected by MIRVA filtering, the mean fre- 
quency is shifted toward the phase frequency. 



IV. APPLICATION TO VORTEX 
FLOWMETERING 

A. Background 

Next, an application of MIRVA filters to vortex 
flowmetering is discussed. Vortex flowmeters are widely 
used in the industry to measure pipe flow. The mea- 
surement principle makes use of the phenomenon of the 
von-Karman vortex street. Behind a shedder bar inserted 
normal to the flow in a pipe, a regular chain of vortices 
is formed, rotating alternatingly left and right. The vol- 
ume flow through the tube can be determined from the 
frequency of vortex formation. In the device used here, 
a piezoelectric sensor sensitive to transversal flow, which 
is inserted downstream behind the shedder bar, is used 
to detect the vortices passing by. A common problem of 
vortex flowmetering is mode locking of the vortex oscilla- 
tions to pulsations in the flow. The second-order statis- 
tics (power spectra) of the sensor signal and the bias on 
the flow measurement in the presence of mode locking 
have been thoroughly investigated [2|j . But it seem pos- 
sible that, by processing the sensor signal with a stronger 
focus on the nonlinear dynamics of the system, a better 
control of mode-locking can be achieved. 

Here we describe the analysis of time series recorded 
in an experiment designed to simulate the problem of de- 
tecting mode locking in an industrial context, using only 
the sensor signal. Methods that have been proposed to 
detect mode locking from univariate time series are the 
analysis of the map of subsequent period lengths of the 
oscillation (angle-of-return-times-map) [3(| and the ap- 
plication of the established bivariate methods on pairs of 
time series extracted from the univariate series by filter- 
ing [3l|. We go along the lines of the second approach, 
making it more powerful by applying MIRVA filters to 
separate the signals. 

The setup of the experiment is sketched in Fig. [3J Pul- 




FIG. 3: Schematic of the experimental setup to record the 
sensor signal of a vortex flowmeter in pulsatile flow. The 
shedder bar has a triangular cross section for efficient vortex 
generation. 
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FIG. 6: The trajectories of demodulated, MIRVA filtered sig- 
nals Zi (see Appendix IA 21 . obtained from Series A at the 
indicated frequencies. The corresponding values of g m i n are 
0.41 0.81 (fpuis), and 0.13 (iw). 



B. Series A: hard lock-in 



FIG. 4: A representative section (top) and the power spec- 
tral density (bottom) of Signal A, which was recorded in the 
experiment sketched in Fig. [3] 
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FIG. 5: The impulse response f(t) of the MIRVA filter cal- 
culated for Series A at the frequency v VO rt- Solid: Re{/(t)}; 
dashed: phase of f(t), relative to an oscillation at constant 
frequency ljo/2-k — 109.64 Hz, i.e., arg [f(t) exp(— iwot)]. The 
overall offset in time is an accidental choice of the search al- 
gorithm. 



sations of the pipe flow were generated by a rotating 
cylinder with three bores orthogonal to the cylinder axis, 
which is inserted into the pipe in such a way that, by 
the rotation, the flow is periodically blocked. This pul- 
sator is driven by an electric motor. The sensor signal 
of a commercial flowmeter, which was mounted about 
40 pipe diameters downstream from the pulsator, was 
recorded. Estimates of the pulsation rate t'puis and the 
frequency of vortex formation f V ort were available on- 
site, while recording the time series. Reynolds numbers 
were O(10 5 ) and the flow was highly turbulent. As a 
result, both inherent and measurement noise contribute 
substantially to the signal. Details of the experimental 
setup will be reported elsewhere. Below, two experimen- 
tal time series labeled as A and B are discussed. Both 
were recorded from the sensor of the vortex flowmeter at 
2 kHz over 250 s. 



1. Extraction of phases 

When recording series A, the flow rate was adjusted 
such as to obtain tVort ~ 110 Hz and the pulsation fre- 
quency was set to f pu i s « f^vort- The time series was 
analysed to determine the strength of the expected 2 : 3 
lock- in. In Fig.0]a representative section of Series A and 
the power spectrum are shown. The oscillations at v VOT t 
can clearly be seen. Since the pulsations themselves do 
not produce any transversal flow, there is only a weak 
signal at v pn i s , presumably due to slight asymmetries in 
the setup. By the nonlinear interaction of vortex street 
and pulsation, flow oscillations at V- = v vmi — v vn \ a 
are excited. These contain significant transversal com- 
ponents and can clearly be seen in the power spectrum. 
The power spectrum also reveals several other oscillatory 
components in the signal. Some of these are nonlinearly 
generated, others are of unknown origin. 

The impulse response of the MIRVA filter for the 
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FIG. 7: Phases, extracted from Series A. Large upper panel: 
the unwrapped phases 3<p-(t) — toot (solid) and <^ V ort(i) — ^ot 
(dashed), where luq/2-k — 109.654Hz (nominal value); 02,3 = 
30- — V ort (dotted). Large lower panel: the cyclic relative 
phase ^2,3- Small panels: Empirical distributions of and 
#2,3 
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110 Hz vortex oscillations (f V ort) is shown in Fig. It 
was calculated by the indirect method described in Ap- 
pendix El using a down-sampling factor h = 50, M = 30 
sampling points, a demodulation frequency luq/2tt — 
109.643 Hz, and a regularization by m = 10th order poly- 
nomials. The overall Gaussian shape of the impulse re- 
sponse of the MIRVA filter can clearly be seen. But 
the filter has additional structure. The oscillation fre- 
quency of the response function decreases with time (see 
the phase in Fig.JHJ. The reason for this particular phase 
dynamics is not clear by now. As can be seen from the 
trajectory of the demodulated filtered signal Zi shown 
in Fig. (right), the phase of the vortex oscillations is 
always well defined. 

From the construction of the pulsator it is clear that 
the flow pulsations have a well defined phase: Each pas- 
sage of a bore along the pulsator inlet (or outlet) defines 
one puis. But the signal-to-noise ratio of the oscillations 
at ^puis is too low to derive unequivocal phase informa- 
tion. As is shown in Fig. [5] (center), the MIRVA filtered 
signal at is pn i s repeatedly reaches the origin of the com- 
plex plane. 

In contrast, the phase of the oscillations at V- is much 
better defined (see Fig. |HJ left) . Since the signal-to-noise 
ratio is smaller at v_ than at v vmt , the MIRVA filter at 
V- is about 8 times more narrow in Fourier space than 
the MIRVA filter at v VOIt . Use of the MIRVA filter (or 
some approximation) is critical for phase extraction at 
2/_ . Here, straightforward boxcar filtering of a region in 
Fourier space containing the V- peak (see, e.g., Ref. |3lj ) 
would be insufficient. 

The phase 0_ of the oscillations at v_ can be used to 
determine the phase pu i s of the pulsator. By the physi- 
cal interpretation of the oscillation at v_ as a nonlinear 
excitation, one has the relation 

4>- = 0vort - ^puls, ( 19 ) 

that yields pu is for known 0_ and </> VO rt- 

2. Relative phases and symmetry 

From the 2 : 3 mode locking, one expects that the 
relative phase 02,3, given by 

4>n,m ■= H^vort ~ ™0puls, (20) 

changes only little over time. For hard mode locking, it 
fluctuates around a constant value. With both hard and 
soft lock-in, the cyclic relative phase ^2.3 defined by 

*n,m := 0n,m mod 2ll (21) 

has an uneven distribution. Generally, one expects the 
distribution to be increasingly sharper localized to a sin- 
gle value as mode locking becomes stronger. The syn- 
chronization index defined in Ref. [12 as 

7^ m : = (cos *„, m } 2 + (sin * n>m > 2 , (22) 



with expectation values estimated by temporal averag- 
ing, was therefore proposed as a measure for the strength 
of mode locking or, more general, phase locking. Absence 
of mode locking is indicated by 7„ jrra = 0, hard coupling 
by 7n,m = 1. In our experiment, we encounter the partic- 
ular situation that vortices and pulsation have opposite 
symmetries with respect to transversal reflection. Thus, 
dynamics is equivariant already with respect to a shift of 
0vort by 7r, rather than 2-7T. Ideally, one would therefore 
always expect 7 f 3 = 0, with or without mode locking. In 
order to take this degeneracy into account, 74^ should be 
used instead of 7 | 3 . 

3. Interpretation of extracted phases 

The evolution of the measured values for 0_ and V ort , 
and of the relative phase 02,3 = 2 V ort — 3 pu is = 
3 4>- ~ 0vort is shown in Fig. [7] Since the definition of 
MIRVA filters leaves the overall delay of the filtered sig- 
nal undetermined, the relative delay of 0_ and V ort has 
to be adjusted in a reasonable way. We choose the de- 
lay such that 72,3 becomes maximal (72,3 turns out not 
to vanish, see below). From the evolution of 02,3 it ap- 
pears that the vortex oscillations contain only a single 
phase slip at about 150 seconds into the time series. But 
upon closer inspection, it appears more plausible to ac- 
count the phase slip to an error in measuring 0_ : As 
expected for this case, the difference in 2 ,3 before and 
after the slip is to a good accuracy 6n. For slips in vor t, 
any other multiple of n would have been possible as well. 
Furthermore, the slip occurs just at the moment when 
the demodulated MIRVA signal at (Fig. El left) goes 
through the small loop reaching toward the coordinate 
center. It appears that the MIRVA filter is too narrow 
for the comparatively low pulsation frequency at this mo- 
ment. In fact, the phase slip disappears when wider filters 
are used - at the price of obtaining new artificial phase 
slips at other times. In conclusion, the data indicate that 
there is not a single real phase slip. Lock-in is hard over 
the full 250 s sampling time. 

The relevant synchronization index 7 | 6 = 0.16 is much 
smaller than one would expect for hard mode locking (see 
also the distributions of ^4,6 in Fig. 0). Use of a syn- 
chronization index based on Shannon entropy [3j yields 
a similar result. Even when the transversal reflection 
symmetry was strongly broken, the then relevant syn- 
chronization index 7 | 3 = 0.62 would be rather low. But, 
as can be seen from Fig. 1101 below, the symmetry is only 
weakly broken. A natural explanation for the discrep- 
ancy between the synchronization index and the phase- 
slip statistics is to assume that most of the fluctuations in 
*2,3, respectively * 4 ,6 (Fig.0 ^4,6 = 2*2,3 mod 2n), are 
due to measurement noise, and not intrinsic to the vortex 
dynamics. This view is compatible with the upper bound 
derived for the variance due to noise in Sec. Ill Bl From 
q- ■= <7max = 0.41 at v- and q VOIt := q max = 0.13 at v VOIt , 
one gets the approximate upper bound 3 g_/4+q£ ort /4 = 
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FIG. 8: The trajectories of demodulated, MIRVA filtered sig- 
nals Zi (see Sec. lA 21 . obtained from Series B and the quotient 
signal x'(t), defined in Sec. II V Cl at the indicated frequencies. 
The corresponding values of q m i n are 0.30 (— ^ pu i s in x'(i)), 
0.39 0_), and 0.22 (iw). 

0.38 for the variance contributed to 4 , 2,3 by measurement 
noise, while the total variance is var\I>2.3 = 0.48. It ap- 
pears that, in some situations, a characterization of mode 
locking by phase-slip statistics is less susceptible to mea- 
surement noise than characterizations based only upon 
cyclic phases, i.e., quantities such as <^ pu is/vort mod 2it or 
- computable thereof - ^4,6- 



C. Series B: soft lock-in 

1. Extraction of phases 

Series B was recorded with iy vor t ~ 110 Hz and ^ pu is ~ 
55 Hz and 1 : 2 mode locking is expected. Similar as 
for Series A, MIRVA filtering readily extracts the un- 
wrapped phase </> V ort of the oscillation ~ exp(i</> vort ) at 
^vort (Fig- HI right). We assume that, as for Series A, 
the oscillations ~ exp(z0 pu i s ) due to pulsation alone are 
much weaker than the oscillations ~ exp(z0 vort — i(j> pu is) 
excited by nonlinear interaction of vortices and pulsation. 
Thus, the nonlinear excitation dominates the oscillations 
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FIG. 9: Unwrapped phases, extracted from Series B: 
20 P uis(t) — coot (solid) and 4> VO it(t) — u;ot (dashed), where 
a>o/27r = 110.697Hz (nominal value); the relative phase 
0i : 2 = 0vort — 20p U i s (dotted). 
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FIG. 10: Empirical distribution function of the cyclic relative 
phase $1,2 obtained from Series B. The double peak is a result 
of the transversal reflection symmetry of the experimental 
system (Fig. 0. 



at v pu \ s w V- := LVort — fpuls- I n contrast to Series A, 
MIRVA filtering at f_ does not lead to an unequivocal 
phase (Fig. [SJ center). It appears that this is due to 
phase slips in <^ vor t, which broaden the range of relevant 
frequencies at z/_ and, as a result, worsen the signal-to- 
noise ratio. 

By making use of the MIRVA filtered signal z V ort(i) 
of the vortex oscillations, </> pu is can nevertheless be ex- 
tracted from the signal. In variation of a method pro- 
posed in Ref. 33], a new, complex- valued time series 
x'(t) := x(t) / z vol: t(t) is constructed from the original 
signal x(t). The overall effect of this transformation is 
to shift all oscillations by f VO rt to negative frequencies. 
The oscillations that were at V- are now at V- — v VOT t = 
— t'puis- They are of the form ~ exp(— z0 pu i s ), i.e., they do 
not depend on the phase of the vortices. Now, MIRVA 
filtering x'(t) at — v pu i s yields the desired unequivocal 
phase information (Fig. [HJ left). As is shown in Fig. EI 
the (fiyort follows <jf> pu is, but several phase slips occur. 

The histogram of the cyclic relative phase ^1,2 reveals 
the two preferred phase angles which are due to transver- 
sal reflection symmetry. Since the symmetry is weakly 
broken, their separation is not exactly tt. Again, it is not 
clear to what degree the broadening of the distribution 
of ^12 is due to measurement noise and to what degree 
to internal noise. 



2. Quantification of the degree of mode locking 



In order to quantify the degree of mode locking inde- 
pendent of measurement noise, a characterization that 
depends only on the long-term dynamics of the phases 
would be useful. Such a measure is, for example, given 
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by ,01,2, with 

_ fl 2 D [(ftyort] + m 2 D [0 pu ls] ~ D [<ftn,m] 

2m 2 D [^puis] 

= lim 

U COV [<ft VQr t(£ + T) - vor t(t), ^ pu l s (£ + T) ~ ^puls(^)] 

2 m var [0 pu i s (* + T) - pu i 8 (t)] 

(23) 

where £)[•] stands for the diffusion coefficient of the spec- 
ified phase variable. /0 n , TO measures in how far the re- 
sponse oscillator (here vortices) follows phase fluctua- 
tions of the drive oscillator (pulsation) . When </> VO rt and 
0puis evolve independently, p n .m = for all n, m. In 
the case of hard n : m lock-in, as was found for Se- 
ries A, D{<j> n ,m] = 0, n 2 D [0vort] = m 2 D [</> pu i s ], and 
Pn,m = 1- Weak mode coupling interpolates between 
these two extremes. Values of p n ,m outside the range 
[0, 1] are possible in principle but unphysical in the sit- 
uation of direct, unidirectional coupling. Since only the 
long-term dynamics of the phases is taken into account, 
rather long time series are required to obtain reproducible 
values of p n , m - For Series B we find, using the estimator 
given by Eq. fBl} with r = 12.5 s, D [0 1)2 ] = 0.4(1) s" 1 , 
Divert] = 2.4(6) s- 1 , and 4£>[0 puIs ] = 2.5(6) s~\ re- 
sulting in pi t 2 ~ 0.9. The empirical value of pi t 2 is stable 
over a wide range in r. Of course, p n ,m could not be used 
when the frequency of the drive oscillator was perfectly 
stable, i.e., when D [4> pu \s] = 0. A detailed analysis of the 
measure p n ,m an d its interpretation is yet to be worked 
out. 



V. CONCLUSION 

MIRVA filtering was introduced as a new method for 
extracting the phase of oscillations from noisy time se- 
ries. It was argued that the phase so obtained is, among 
other favorable properties, particularly robust to noise 
and linear filtering of the signal. Detailed directions for 
computing MIRVA filters numerically are given in Ap- 
pendix^ In a numerical case study it was demonstrated 
that MIRVA filtering introduces only little bias to the 
phase (or average) frequency. A new synchronization in- 
dex has been proposed, which is designed to be robust to 
noise if MIRVA filtering is used. 

By applying MIRVA filtering to the signal of a vortex 
flowmeter, we showed that the method can be used to ob- 
tain well defined phases from oscillatory time series even 
under unfavorable conditions such as strong internal and 
measurement noise. The phases were used to investigate 
the strength of mode locking. As another application, 
the MIRVA filtered signal was used for a data driven de- 
modulation technique in Sec. II V CI Limitations to phase 
extraction, which remain even when MIRVA filters are 
used, have been addressed. 
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APPENDIX A: COMPUTATION OF MIRVA 
FILTERS 

1. Main algorithms 

In practical applications time series Xi (i = 1, ...,JV), 
sampled from x(t) at evenly spaced discrete times t = iAt 
are given. The MIRVA filters have to be estimated from 
this data. Here, two methods are proposed. The first 
method is more appropriate for short time series, the sec- 
ond is more efficient when N is large. With both meth- 
ods, the impulse response fj of the filters is restricted to 
a finite length M (j = 1, . . . , M). 



a. Direct method 
When using the first method, the convolution 

M 

Zk = (x* f)k+M = ^ fj x k+M-j (Al) 
3=1 

(for k = 1, . . . , N — M+ 1) is calculated directly, and the 
expectation values in Eq. are estimated as averages 
over k. An iterative search algorithm is used to find the 
MIRVA filter fj with q 2 = q 2 min . 

b. Indirect method 

The second method makes use of the fact that q de- 
pends on x(t) only through its second and fourth mo- 
ment. It is often more efficient than the first method 
but, as a trade off, entails systematic errors of the order 
0(M/N) in the estimation of the moments of Zi. For 
notational convenience, we define the second and fourth 
moments Xi in the time-discrete representation as 

C-ijkl — (%r—i X r —j X T — k X T —lj (A2) 

and 

C{ j (x T _^ x T —j^ (A3) 

with arbitrary r. These expectation values are estimated 
by averaging over time and making use of symmetries, 
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Equation JQ now reads 

2 Sij fci C »j klfifjfkfl , . ., 

9 = —7 —2 1 ( A4 ) 

\^>ij °i jf if j J 

with all sums running over 1, . . . , M. Thus, while 
the computation of the moments takes time of order 
0(N M 3 ), the time required for the optimization itself 
is independent of N. Besides, derivatives of q 2 with re- 
spect to ft are calculated at little additional cost, and 
can be provided to the optimization algorithm to help 
finding a MIRVA filter f t . 



2. Down-sampling, demodulation, and 
regularization 

In order not to introduce artificial restrictions of the 
search space for /j, the duration of the impulse response, 
i.e. M x At should be of the order of the phase coherence 
time of the oscillation. This time can easily cover several 
hundred oscillation periods. Computation time depends 
critically on M. To keep M low and make the compu- 
tation feasible, it is therefore advisable not to work with 
the raw time series Xi but to perform a demodulation 
and a down-sampling step prior to the main calculation. 
For the calculations discussed in Sees. IIIIIIVI instead of 
Xi, the demodulated time series Xj given by 

Xj := ^2 K i x hj+i exp[-i (hj + I) uj At], (A5) 
l 

were used with a symmetric, triangular smoothing kernel 
Ki at a width of two time the down-sampling factor h. 
The demodulation frequency luq should be set to a value 
close to frequency of the targeted oscillations. 

To see the effect of this transformation, notice that for 
stationary, discrete-time processes the value of q defined 
by Eq. (|TJ with z replaced by Z — F * X is for any filter 
Fk identical to the value obtained with z — f*x, provided 

f k = ^Ki F {k+l)/h exp(i k w At) (A6) 
i 

and Fk := for non-integer k by convention: It is not 
difficult to verify that Zj = Zhj exp(— ihjujoAt), indepen- 
dent of K[. 

As a result, every MIRVA filter Fk for Xj leads by 
Eq. (|A6|I to the approximate MIRVA filter f k for Xi . The 
approximation is good if the interpolation 1A6(I of Fk 
defined by Ki is reasonable. 

The MIRVA filters Fk found for typical experimental 
data are more or less deformed variants of Gaussian filters 
F k w exp(-±(fc-M/2) 2 /i 2 Aw 2 At 2 ) with bandwidth Alu. 
The linear interpolation for fk given by the triangular K\ 
is good if h Auj At <C 1 . In practice, this requires filter 
lengths of at least M * 15 — 30. 



Experimental time series are often not long enough to 
yield faithful estimates for all 0(M 3 ) independent ele- 
ments of Cijki- This problem can be overcome by a 
regularization of Fk- In our calculations, we restricted 
the filters Fk to the family Fk — exp(P TO (fc)), with m-th 
order polynomials P m . 



3. A guide to choosing appropriate parameters 

The following procedures were used to find appropriate 
values for the demodulation frequency luq and the down- 
sampling factor h, which determines the duration of the 
filter's impulse response hMAt. The sampling rate At is 
assumed to be given and M is restricted to a small range 
by computational limitations. 

After an initial guess, the value of u>o was set to the 
value of the empirically found phase frequency uj p h of Zi 
in an iterative process. In order to adjust h, the enve- 
lope \Fk \ of the computed MIRVA filter was investigated. 
When h is too large, \Fk\ has a sharp peak and vanishes 
for all other values. When h is too small, most weight of 
the filter is concentrated near the end points F% and Fm- 
By inspection one finds F\ s» ±iFM, i.e., the MIRVA 
filter with a constraint in the filter length approximates 
a simple 2D delay embedding, h is adjusted accordingly. 

A systematic procedure for finding good values for the 
polynomial order m has not been developed yet. But 
with Tn ~ 6 — 10 results generally depend little on the 
precise value. 



4. Convergence and side minima 

In Sec. ^ it was proposed to identify local minima of 
q with distinct oscillatory components of the signal. The 
structure of the search space is therefore of immediate 
theoretical interest. In fact, with long enough time series 
sampled from a typical signal, the numerical search algo- 
rithms used here consistently and effortlessly reach local 
minima located in a small set of well separated points in 
the space of all filters, irrespective of the - randomly cho- 
sen - starting points. A unique minimum is typically sin- 
gled out when using demodulation and down sampling, 
since this effectively implies a pre-selection of the fre- 
quency range of interest. With shorter time series, how- 
ever, these minima split into large clusters of several side 
minima. In order to cope with these artificial multiplici- 
ties, only those local minima were accepted for selecting 
MIRVA filters which where found three times within a 
series of minimization runs with random starting points, 
without previously finding any point with a lower value 
of q. 
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APPENDIX B: REMARKS ON THE 
ESTIMATION OF D 

We discuss a method to determine the diffusion coef- 
ficient D defined by Eq. i|14fl from samples 4>(t) of fi- 
nite length (0 < t < T). First w p h is estimated as 
i2) p h = [(A(T) - (f)(0)] /T. Then D can be estimated by 



D T := 



/ T - T [0( t + T)-#t)-a 



phTj 



2t(T-t)(1 -r/T) 



(Bl) 



where < r < T. The last factor in the denominator 
compensates for the loss of statistical degrees of freedom 
by the estimation of o> p h as Co v h_. When assuming 4>(t) to 
perform a random walk with constant drift, it is straight- 
forwardly verified that tD p h is a maximum likelihood es- 
timator and (D r ) = D. Under the same assumption, the 
variance of u) p h is 2D/T and 

cov(D k T,DiT) = 



D 



Z[6(l - kfk - 2 (1 - k) (1 + 3/c 2 ) I - (1 - Ak) I 2 



3fc(l-fc) 2 (l-0 2 



(B2) 



for 0<fc<Z<l/2 (the last expression was obtained 
with the help of symbolic computer algebra). In partic- 
ular, 



var Dit = D 



I (4- lll + Al 2 



6 1 3 



3(1-0' 



(B3) 



which increases monotonically with < I < 1/2. For 
good estimates of D, the parameter I should be chosen 
as small as possible but large enough to justify the as- 
sumption of a random walk over times I T . The estimator 
D T can be slightly improved by using linear combinations 
with different r. For example, the variance of 



DL 



-Dr 



(B4) 



is about 10% smaller than of D T , as is verified using 
Eq. CE2J. 
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